using System;

namespace ICP_MS.Plugin.HPLC;

public class LEquations
{
	private Matrix mtxLECoef;

	private Matrix mtxLEConst;

	public LEquations()
	{
	}

	public LEquations(Matrix mtxCoef, Matrix mtxConst)
	{
		Init(mtxCoef, mtxConst);
	}

	public bool Init(Matrix mtxCoef, Matrix mtxConst)
	{
		if (mtxCoef.GetNumRows() != mtxConst.GetNumRows())
		{
			return false;
		}
		mtxLECoef = new Matrix(mtxCoef);
		mtxLEConst = new Matrix(mtxConst);
		return true;
	}

	public Matrix GetCoefMatrix()
	{
		return mtxLECoef;
	}

	public Matrix GetConstMatrix()
	{
		return mtxLEConst;
	}

	public int GetNumEquations()
	{
		return GetCoefMatrix().GetNumRows();
	}

	public int GetNumUnknowns()
	{
		return GetCoefMatrix().GetNumColumns();
	}

	public bool GetRootsetGauss(Matrix mtxResult)
	{
		int nIs = 0;
		mtxResult.SetValue(mtxLEConst);
		double[] pDataCoef = mtxLECoef.GetData();
		double[] pDataConst = mtxResult.GetData();
		int m = GetNumUnknowns();
		int[] pnJs = new int[m];
		int l = 1;
		double d;
		for (int k = 0; k <= m - 2; k++)
		{
			d = 0.0;
			for (int i = k; i <= m - 1; i++)
			{
				for (int j = k; j <= m - 1; j++)
				{
					double t = Math.Abs(pDataCoef[i * m + j]);
					if (t > d)
					{
						d = t;
						pnJs[k] = j;
						nIs = i;
					}
				}
			}
			if (d == 0.0)
			{
				l = 0;
			}
			else
			{
				if (pnJs[k] != k)
				{
					for (int i = 0; i <= m - 1; i++)
					{
						int p = i * m + k;
						int q = i * m + pnJs[k];
						double t = pDataCoef[p];
						pDataCoef[p] = pDataCoef[q];
						pDataCoef[q] = t;
					}
				}
				if (nIs != k)
				{
					double t;
					for (int j = k; j <= m - 1; j++)
					{
						int p = k * m + j;
						int q = nIs * m + j;
						t = pDataCoef[p];
						pDataCoef[p] = pDataCoef[q];
						pDataCoef[q] = t;
					}
					t = pDataConst[k];
					pDataConst[k] = pDataConst[nIs];
					pDataConst[nIs] = t;
				}
			}
			if (l == 0)
			{
				return false;
			}
			d = pDataCoef[k * m + k];
			for (int j = k + 1; j <= m - 1; j++)
			{
				int p = k * m + j;
				pDataCoef[p] /= d;
			}
			pDataConst[k] /= d;
			for (int i = k + 1; i <= m - 1; i++)
			{
				for (int j = k + 1; j <= m - 1; j++)
				{
					int p = i * m + j;
					pDataCoef[p] -= pDataCoef[i * m + k] * pDataCoef[k * m + j];
				}
				pDataConst[i] -= pDataCoef[i * m + k] * pDataConst[k];
			}
		}
		d = pDataCoef[(m - 1) * m + m - 1];
		if (d == 0.0)
		{
			return false;
		}
		pDataConst[m - 1] /= d;
		for (int i = m - 2; i >= 0; i--)
		{
			double t = 0.0;
			for (int j = i + 1; j <= m - 1; j++)
			{
				t += pDataCoef[i * m + j] * pDataConst[j];
			}
			pDataConst[i] -= t;
		}
		pnJs[m - 1] = m - 1;
		for (int k = m - 1; k >= 0; k--)
		{
			if (pnJs[k] != k)
			{
				double t = pDataConst[k];
				pDataConst[k] = pDataConst[pnJs[k]];
				pDataConst[pnJs[k]] = t;
			}
		}
		return true;
	}

	public bool GetRootsetGaussJordan(Matrix mtxResult)
	{
		int nIs = 0;
		mtxResult.SetValue(mtxLEConst);
		double[] pDataCoef = mtxLECoef.GetData();
		double[] pDataConst = mtxResult.GetData();
		int n = GetNumUnknowns();
		int m = mtxLEConst.GetNumColumns();
		int[] pnJs = new int[n];
		int l = 1;
		for (int k = 0; k <= n - 1; k++)
		{
			double d = 0.0;
			for (int i = k; i <= n - 1; i++)
			{
				for (int j = k; j <= n - 1; j++)
				{
					double t = Math.Abs(pDataCoef[i * n + j]);
					if (t > d)
					{
						d = t;
						pnJs[k] = j;
						nIs = i;
					}
				}
			}
			if (d + 1.0 == 1.0)
			{
				l = 0;
			}
			else
			{
				if (pnJs[k] != k)
				{
					for (int i = 0; i <= n - 1; i++)
					{
						int p = i * n + k;
						int q = i * n + pnJs[k];
						double t = pDataCoef[p];
						pDataCoef[p] = pDataCoef[q];
						pDataCoef[q] = t;
					}
				}
				if (nIs != k)
				{
					for (int j = k; j <= n - 1; j++)
					{
						int p = k * n + j;
						int q = nIs * n + j;
						double t = pDataCoef[p];
						pDataCoef[p] = pDataCoef[q];
						pDataCoef[q] = t;
					}
					for (int j = 0; j <= m - 1; j++)
					{
						int p = k * m + j;
						int q = nIs * m + j;
						double t = pDataConst[p];
						pDataConst[p] = pDataConst[q];
						pDataConst[q] = t;
					}
				}
			}
			if (l == 0)
			{
				return false;
			}
			d = pDataCoef[k * n + k];
			for (int j = k + 1; j <= n - 1; j++)
			{
				int p = k * n + j;
				pDataCoef[p] /= d;
			}
			for (int j = 0; j <= m - 1; j++)
			{
				int p = k * m + j;
				pDataConst[p] /= d;
			}
			for (int j = k + 1; j <= n - 1; j++)
			{
				for (int i = 0; i <= n - 1; i++)
				{
					int p = i * n + j;
					if (i != k)
					{
						pDataCoef[p] -= pDataCoef[i * n + k] * pDataCoef[k * n + j];
					}
				}
			}
			for (int j = 0; j <= m - 1; j++)
			{
				for (int i = 0; i <= n - 1; i++)
				{
					int p = i * m + j;
					if (i != k)
					{
						pDataConst[p] -= pDataCoef[i * n + k] * pDataConst[k * m + j];
					}
				}
			}
		}
		for (int k = n - 1; k >= 0; k--)
		{
			if (pnJs[k] != k)
			{
				for (int j = 0; j <= m - 1; j++)
				{
					int p = k * m + j;
					int q = pnJs[k] * m + j;
					double t = pDataConst[p];
					pDataConst[p] = pDataConst[q];
					pDataConst[q] = t;
				}
			}
		}
		return true;
	}

	public bool GetRootsetGauss(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag)
	{
		int nIs = 0;
		mtxResult.SetValue(mtxLEConst);
		mtxResultImag.SetValue(mtxConstImag);
		double[] pDataCoef = mtxLECoef.GetData();
		double[] pDataConst = mtxResult.GetData();
		double[] pDataCoefImag = mtxCoefImag.GetData();
		double[] pDataConstImag = mtxResultImag.GetData();
		int n = GetNumUnknowns();
		int m = mtxLEConst.GetNumColumns();
		int[] pnJs = new int[n];
		double d;
		int u;
		double q;
		double s;
		double p;
		for (int k = 0; k <= n - 2; k++)
		{
			d = 0.0;
			for (int i = k; i <= n - 1; i++)
			{
				for (int j = k; j <= n - 1; j++)
				{
					u = i * n + j;
					p = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];
					if (p > d)
					{
						d = p;
						pnJs[k] = j;
						nIs = i;
					}
				}
			}
			if (d == 0.0)
			{
				return false;
			}
			int v;
			if (nIs != k)
			{
				for (int j = k; j <= n - 1; j++)
				{
					u = k * n + j;
					v = nIs * n + j;
					p = pDataCoef[u];
					pDataCoef[u] = pDataCoef[v];
					pDataCoef[v] = p;
					p = pDataCoefImag[u];
					pDataCoefImag[u] = pDataCoefImag[v];
					pDataCoefImag[v] = p;
				}
				p = pDataConst[k];
				pDataConst[k] = pDataConst[nIs];
				pDataConst[nIs] = p;
				p = pDataConstImag[k];
				pDataConstImag[k] = pDataConstImag[nIs];
				pDataConstImag[nIs] = p;
			}
			if (pnJs[k] != k)
			{
				for (int i = 0; i <= n - 1; i++)
				{
					u = i * n + k;
					v = i * n + pnJs[k];
					p = pDataCoef[u];
					pDataCoef[u] = pDataCoef[v];
					pDataCoef[v] = p;
					p = pDataCoefImag[u];
					pDataCoefImag[u] = pDataCoefImag[v];
					pDataCoefImag[v] = p;
				}
			}
			v = k * n + k;
			for (int j = k + 1; j <= n - 1; j++)
			{
				u = k * n + j;
				p = pDataCoef[u] * pDataCoef[v];
				q = (0.0 - pDataCoefImag[u]) * pDataCoefImag[v];
				s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataCoef[u] + pDataCoefImag[u]);
				pDataCoef[u] = (p - q) / d;
				pDataCoefImag[u] = (s - p - q) / d;
			}
			p = pDataConst[k] * pDataCoef[v];
			q = (0.0 - pDataConstImag[k]) * pDataCoefImag[v];
			s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataConst[k] + pDataConstImag[k]);
			pDataConst[k] = (p - q) / d;
			pDataConstImag[k] = (s - p - q) / d;
			for (int i = k + 1; i <= n - 1; i++)
			{
				u = i * n + k;
				for (int j = k + 1; j <= n - 1; j++)
				{
					v = k * n + j;
					int l = i * n + j;
					p = pDataCoef[u] * pDataCoef[v];
					q = pDataCoefImag[u] * pDataCoefImag[v];
					s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataCoef[v] + pDataCoefImag[v]);
					pDataCoef[l] = pDataCoef[l] - p + q;
					pDataCoefImag[l] = pDataCoefImag[l] - s + p + q;
				}
				p = pDataCoef[u] * pDataConst[k];
				q = pDataCoefImag[u] * pDataConstImag[k];
				s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[k] + pDataConstImag[k]);
				pDataConst[i] = pDataConst[i] - p + q;
				pDataConstImag[i] = pDataConstImag[i] - s + p + q;
			}
		}
		u = (n - 1) * n + n - 1;
		d = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];
		if (d == 0.0)
		{
			return false;
		}
		p = pDataCoef[u] * pDataConst[n - 1];
		q = (0.0 - pDataCoefImag[u]) * pDataConstImag[n - 1];
		s = (pDataCoef[u] - pDataCoefImag[u]) * (pDataConst[n - 1] + pDataConstImag[n - 1]);
		pDataConst[n - 1] = (p - q) / d;
		pDataConstImag[n - 1] = (s - p - q) / d;
		for (int i = n - 2; i >= 0; i--)
		{
			for (int j = i + 1; j <= n - 1; j++)
			{
				u = i * n + j;
				p = pDataCoef[u] * pDataConst[j];
				q = pDataCoefImag[u] * pDataConstImag[j];
				s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[j] + pDataConstImag[j]);
				pDataConst[i] = pDataConst[i] - p + q;
				pDataConstImag[i] = pDataConstImag[i] - s + p + q;
			}
		}
		pnJs[n - 1] = n - 1;
		for (int k = n - 1; k >= 0; k--)
		{
			if (pnJs[k] != k)
			{
				p = pDataConst[k];
				pDataConst[k] = pDataConst[pnJs[k]];
				pDataConst[pnJs[k]] = p;
				p = pDataConstImag[k];
				pDataConstImag[k] = pDataConstImag[pnJs[k]];
				pDataConstImag[pnJs[k]] = p;
			}
		}
		return true;
	}

	public bool GetRootsetGaussJordan(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag)
	{
		int nIs = 0;
		mtxResult.SetValue(mtxLEConst);
		mtxResultImag.SetValue(mtxConstImag);
		double[] pDataCoef = mtxLECoef.GetData();
		double[] pDataConst = mtxResult.GetData();
		double[] pDataCoefImag = mtxCoefImag.GetData();
		double[] pDataConstImag = mtxResultImag.GetData();
		int n = GetNumUnknowns();
		int m = mtxLEConst.GetNumColumns();
		int[] pnJs = new int[n];
		for (int k = 0; k <= n - 1; k++)
		{
			double d = 0.0;
			for (int i = k; i <= n - 1; i++)
			{
				for (int j = k; j <= n - 1; j++)
				{
					int u = i * n + j;
					double p = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];
					if (p > d)
					{
						d = p;
						pnJs[k] = j;
						nIs = i;
					}
				}
			}
			if (d == 0.0)
			{
				return false;
			}
			int v;
			if (nIs != k)
			{
				for (int j = k; j <= n - 1; j++)
				{
					int u = k * n + j;
					v = nIs * n + j;
					double p = pDataCoef[u];
					pDataCoef[u] = pDataCoef[v];
					pDataCoef[v] = p;
					p = pDataCoefImag[u];
					pDataCoefImag[u] = pDataCoefImag[v];
					pDataCoefImag[v] = p;
				}
				for (int j = 0; j <= m - 1; j++)
				{
					int u = k * m + j;
					v = nIs * m + j;
					double p = pDataConst[u];
					pDataConst[u] = pDataConst[v];
					pDataConst[v] = p;
					p = pDataConstImag[u];
					pDataConstImag[u] = pDataConstImag[v];
					pDataConstImag[v] = p;
				}
			}
			if (pnJs[k] != k)
			{
				for (int i = 0; i <= n - 1; i++)
				{
					int u = i * n + k;
					v = i * n + pnJs[k];
					double p = pDataCoef[u];
					pDataCoef[u] = pDataCoef[v];
					pDataCoef[v] = p;
					p = pDataCoefImag[u];
					pDataCoefImag[u] = pDataCoefImag[v];
					pDataCoefImag[v] = p;
				}
			}
			v = k * n + k;
			for (int j = k + 1; j <= n - 1; j++)
			{
				int u = k * n + j;
				double p = pDataCoef[u] * pDataCoef[v];
				double q = (0.0 - pDataCoefImag[u]) * pDataCoefImag[v];
				double s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataCoef[u] + pDataCoefImag[u]);
				pDataCoef[u] = (p - q) / d;
				pDataCoefImag[u] = (s - p - q) / d;
			}
			for (int j = 0; j <= m - 1; j++)
			{
				int u = k * m + j;
				double p = pDataConst[u] * pDataCoef[v];
				double q = (0.0 - pDataConstImag[u]) * pDataCoefImag[v];
				double s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataConst[u] + pDataConstImag[u]);
				pDataConst[u] = (p - q) / d;
				pDataConstImag[u] = (s - p - q) / d;
			}
			for (int i = 0; i <= n - 1; i++)
			{
				if (i != k)
				{
					int u = i * n + k;
					for (int j = k + 1; j <= n - 1; j++)
					{
						v = k * n + j;
						int l = i * n + j;
						double p = pDataCoef[u] * pDataCoef[v];
						double q = pDataCoefImag[u] * pDataCoefImag[v];
						double s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataCoef[v] + pDataCoefImag[v]);
						pDataCoef[l] = pDataCoef[l] - p + q;
						pDataCoefImag[l] = pDataCoefImag[l] - s + p + q;
					}
					for (int j = 0; j <= m - 1; j++)
					{
						int l = i * m + j;
						v = k * m + j;
						double p = pDataCoef[u] * pDataConst[v];
						double q = pDataCoefImag[u] * pDataConstImag[v];
						double s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[v] + pDataConstImag[v]);
						pDataConst[l] = pDataConst[l] - p + q;
						pDataConstImag[l] = pDataConstImag[l] - s + p + q;
					}
				}
			}
		}
		for (int k = n - 1; k >= 0; k--)
		{
			if (pnJs[k] != k)
			{
				for (int j = 0; j <= m - 1; j++)
				{
					int u = k * m + j;
					int v = pnJs[k] * m + j;
					double p = pDataConst[u];
					pDataConst[u] = pDataConst[v];
					pDataConst[v] = p;
					p = pDataConstImag[u];
					pDataConstImag[u] = pDataConstImag[v];
					pDataConstImag[v] = p;
				}
			}
		}
		return true;
	}

	public bool GetRootsetTriDiagonal(Matrix mtxResult)
	{
		mtxResult.SetValue(mtxLEConst);
		double[] pDataConst = mtxResult.GetData();
		int k = GetNumUnknowns();
		if (mtxLECoef.GetNumRows() != k)
		{
			return false;
		}
		double[] pDiagData = new double[3 * k - 2];
		int i;
		int j = (i = 0);
		if (j == 0)
		{
			pDiagData[i++] = mtxLECoef.GetElement(j, j);
			pDiagData[i++] = mtxLECoef.GetElement(j, j + 1);
		}
		for (j = 1; j < k - 1; j++)
		{
			pDiagData[i++] = mtxLECoef.GetElement(j, j - 1);
			pDiagData[i++] = mtxLECoef.GetElement(j, j);
			pDiagData[i++] = mtxLECoef.GetElement(j, j + 1);
		}
		if (j == k - 1)
		{
			pDiagData[i++] = mtxLECoef.GetElement(j, j - 1);
			pDiagData[i++] = mtxLECoef.GetElement(j, j);
		}
		double s;
		for (j = 0; j <= k - 2; j++)
		{
			i = 3 * j;
			s = pDiagData[i];
			if (Math.Abs(s) + 1.0 == 1.0)
			{
				return false;
			}
			pDiagData[i + 1] /= s;
			pDataConst[j] /= s;
			pDiagData[i + 3] -= pDiagData[i + 2] * pDiagData[i + 1];
			pDataConst[j + 1] -= pDiagData[i + 2] * pDataConst[j];
		}
		s = pDiagData[3 * k - 3];
		if (s == 0.0)
		{
			return false;
		}
		pDataConst[k - 1] /= s;
		for (j = k - 2; j >= 0; j--)
		{
			pDataConst[j] -= pDiagData[3 * j + 1] * pDataConst[j + 1];
		}
		return true;
	}

	public bool GetRootsetBand(int nBandWidth, Matrix mtxResult)
	{
		int nis = 0;
		if ((nBandWidth - 1) % 2 != 0)
		{
			return false;
		}
		mtxResult.SetValue(mtxLEConst);
		double[] pDataConst = mtxResult.GetData();
		int l = mtxLEConst.GetNumColumns();
		int m = GetNumUnknowns();
		if (mtxLECoef.GetNumRows() != m)
		{
			return false;
		}
		double[] pBandData = new double[nBandWidth * m];
		int ls = (nBandWidth - 1) / 2;
		for (int i = 0; i < m; i++)
		{
			int j = 0;
			for (int k = Math.Max(0, i - ls); k < Math.Max(0, i - ls) + nBandWidth; k++)
			{
				if (k < m)
				{
					pBandData[i * nBandWidth + j++] = mtxLECoef.GetElement(i, k);
				}
				else
				{
					pBandData[i * nBandWidth + j++] = 0.0;
				}
			}
		}
		double p;
		for (int k = 0; k <= m - 2; k++)
		{
			p = 0.0;
			for (int i = k; i <= ls; i++)
			{
				double t = Math.Abs(pBandData[i * nBandWidth]);
				if (t > p)
				{
					p = t;
					nis = i;
				}
			}
			if (p == 0.0)
			{
				return false;
			}
			for (int j = 0; j <= l - 1; j++)
			{
				int u = k * l + j;
				int v = nis * l + j;
				double t = pDataConst[u];
				pDataConst[u] = pDataConst[v];
				pDataConst[v] = t;
			}
			for (int j = 0; j <= nBandWidth - 1; j++)
			{
				int u = k * nBandWidth + j;
				int v = nis * nBandWidth + j;
				double t = pBandData[u];
				pBandData[u] = pBandData[v];
				pBandData[v] = t;
			}
			for (int j = 0; j <= l - 1; j++)
			{
				int u = k * l + j;
				pDataConst[u] /= pBandData[k * nBandWidth];
			}
			for (int j = 1; j <= nBandWidth - 1; j++)
			{
				int u = k * nBandWidth + j;
				pBandData[u] /= pBandData[k * nBandWidth];
			}
			for (int i = k + 1; i <= ls; i++)
			{
				double t = pBandData[i * nBandWidth];
				int u;
				for (int j = 0; j <= l - 1; j++)
				{
					u = i * l + j;
					int v = k * l + j;
					pDataConst[u] -= t * pDataConst[v];
				}
				for (int j = 1; j <= nBandWidth - 1; j++)
				{
					u = i * nBandWidth + j;
					int v = k * nBandWidth + j;
					pBandData[u - 1] = pBandData[u] - t * pBandData[v];
				}
				u = i * nBandWidth + nBandWidth - 1;
				pBandData[u] = 0.0;
			}
			if (ls != m - 1)
			{
				ls++;
			}
		}
		p = pBandData[(m - 1) * nBandWidth];
		if (p == 0.0)
		{
			return false;
		}
		for (int j = 0; j <= l - 1; j++)
		{
			int u = (m - 1) * l + j;
			pDataConst[u] /= p;
		}
		ls = 1;
		for (int i = m - 2; i >= 0; i--)
		{
			for (int k = 0; k <= l - 1; k++)
			{
				int u = i * l + k;
				for (int j = 1; j <= ls; j++)
				{
					int v = i * nBandWidth + j;
					nis = (i + j) * l + k;
					pDataConst[u] -= pBandData[v] * pDataConst[nis];
				}
			}
			if (ls != nBandWidth - 1)
			{
				ls++;
			}
		}
		return true;
	}

	public bool GetRootsetDjn(Matrix mtxResult)
	{
		Matrix mtxCoef = new Matrix(mtxLECoef);
		mtxResult.SetValue(mtxLEConst);
		int n = mtxCoef.GetNumColumns();
		int m = mtxResult.GetNumColumns();
		double[] pDataCoef = mtxCoef.GetData();
		double[] pDataConst = mtxResult.GetData();
		if (pDataCoef[0] == 0.0)
		{
			return false;
		}
		int u;
		for (int i = 1; i <= n - 1; i++)
		{
			u = i * n;
			pDataCoef[u] /= pDataCoef[0];
		}
		double p;
		for (int i = 1; i <= n - 2; i++)
		{
			u = i * n + i;
			for (int j = 1; j <= i; j++)
			{
				int v = i * n + j - 1;
				int l = (j - 1) * n + j - 1;
				pDataCoef[u] -= pDataCoef[v] * pDataCoef[v] * pDataCoef[l];
			}
			p = pDataCoef[u];
			if (p == 0.0)
			{
				return false;
			}
			for (int k = i + 1; k <= n - 1; k++)
			{
				u = k * n + i;
				for (int j = 1; j <= i; j++)
				{
					int v = k * n + j - 1;
					int l = i * n + j - 1;
					int w = (j - 1) * n + j - 1;
					pDataCoef[u] -= pDataCoef[v] * pDataCoef[l] * pDataCoef[w];
				}
				pDataCoef[u] /= p;
			}
		}
		u = n * n - 1;
		for (int j = 1; j <= n - 1; j++)
		{
			int v = (n - 1) * n + j - 1;
			int w = (j - 1) * n + j - 1;
			pDataCoef[u] -= pDataCoef[v] * pDataCoef[v] * pDataCoef[w];
		}
		p = pDataCoef[u];
		if (p == 0.0)
		{
			return false;
		}
		for (int j = 0; j <= m - 1; j++)
		{
			for (int i = 1; i <= n - 1; i++)
			{
				u = i * m + j;
				for (int k = 1; k <= i; k++)
				{
					int v = i * n + k - 1;
					int w = (k - 1) * m + j;
					pDataConst[u] -= pDataCoef[v] * pDataConst[w];
				}
			}
		}
		for (int i = 1; i <= n - 1; i++)
		{
			u = (i - 1) * n + i - 1;
			for (int j = i; j <= n - 1; j++)
			{
				int v = (i - 1) * n + j;
				int w = j * n + i - 1;
				pDataCoef[v] = pDataCoef[u] * pDataCoef[w];
			}
		}
		for (int j = 0; j <= m - 1; j++)
		{
			u = (n - 1) * m + j;
			pDataConst[u] /= p;
			for (int k = 1; k <= n - 1; k++)
			{
				int k2 = n - k;
				int k4 = k2 - 1;
				u = k4 * m + j;
				for (int k3 = k2; k3 <= n - 1; k3++)
				{
					int v = k4 * n + k3;
					int w = k3 * m + j;
					pDataConst[u] -= pDataCoef[v] * pDataConst[w];
				}
				pDataConst[u] /= pDataCoef[k4 * n + k4];
			}
		}
		return true;
	}

	public bool GetRootsetCholesky(Matrix mtxResult)
	{
		Matrix mtxCoef = new Matrix(mtxLECoef);
		mtxResult.SetValue(mtxLEConst);
		int m = mtxCoef.GetNumColumns();
		int l = mtxResult.GetNumColumns();
		double[] pDataCoef = mtxCoef.GetData();
		double[] pDataConst = mtxResult.GetData();
		if (pDataCoef[0] <= 0.0)
		{
			return false;
		}
		pDataCoef[0] = Math.Sqrt(pDataCoef[0]);
		for (int j = 1; j <= m - 1; j++)
		{
			pDataCoef[j] /= pDataCoef[0];
		}
		for (int i = 1; i <= m - 1; i++)
		{
			int u = i * m + i;
			for (int j = 1; j <= i; j++)
			{
				int v = (j - 1) * m + i;
				pDataCoef[u] -= pDataCoef[v] * pDataCoef[v];
			}
			if (pDataCoef[u] <= 0.0)
			{
				return false;
			}
			pDataCoef[u] = Math.Sqrt(pDataCoef[u]);
			if (i == m - 1)
			{
				continue;
			}
			for (int j = i + 1; j <= m - 1; j++)
			{
				int v = i * m + j;
				for (int k = 1; k <= i; k++)
				{
					pDataCoef[v] -= pDataCoef[(k - 1) * m + i] * pDataCoef[(k - 1) * m + j];
				}
				pDataCoef[v] /= pDataCoef[u];
			}
		}
		for (int j = 0; j <= l - 1; j++)
		{
			pDataConst[j] /= pDataCoef[0];
			for (int i = 1; i <= m - 1; i++)
			{
				int u = i * m + i;
				int v = i * l + j;
				for (int k = 1; k <= i; k++)
				{
					pDataConst[v] -= pDataCoef[(k - 1) * m + i] * pDataConst[(k - 1) * l + j];
				}
				pDataConst[v] /= pDataCoef[u];
			}
		}
		for (int j = 0; j <= l - 1; j++)
		{
			int u = (m - 1) * l + j;
			pDataConst[u] /= pDataCoef[m * m - 1];
			for (int k = m - 1; k >= 1; k--)
			{
				u = (k - 1) * l + j;
				int v;
				for (int i = k; i <= m - 1; i++)
				{
					v = (k - 1) * m + i;
					pDataConst[u] -= pDataCoef[v] * pDataConst[i * l + j];
				}
				v = (k - 1) * m + k - 1;
				pDataConst[u] /= pDataCoef[v];
			}
		}
		return true;
	}

	public bool GetRootsetGgje(Matrix mtxResult)
	{
		int nIs = 0;
		Matrix mtxCoef = new Matrix(mtxLECoef);
		mtxResult.SetValue(mtxLEConst);
		int l = mtxCoef.GetNumColumns();
		double[] pDataCoef = mtxCoef.GetData();
		double[] pDataConst = mtxResult.GetData();
		int[] pnJs = new int[l];
		for (int k = 0; k <= l - 1; k++)
		{
			double d = 0.0;
			double t;
			for (int i = k; i <= l - 1; i++)
			{
				for (int j = k; j <= l - 1; j++)
				{
					t = Math.Abs(pDataCoef[i * l + j]);
					if (t > d)
					{
						d = t;
						pnJs[k] = j;
						nIs = i;
					}
				}
			}
			if (d == 0.0)
			{
				return false;
			}
			if (nIs != k)
			{
				for (int j = k; j <= l - 1; j++)
				{
					int u = k * l + j;
					int v = nIs * l + j;
					t = pDataCoef[u];
					pDataCoef[u] = pDataCoef[v];
					pDataCoef[v] = t;
				}
				t = pDataConst[k];
				pDataConst[k] = pDataConst[nIs];
				pDataConst[nIs] = t;
			}
			if (pnJs[k] != k)
			{
				for (int i = 0; i <= l - 1; i++)
				{
					int u = i * l + k;
					int v = i * l + pnJs[k];
					t = pDataCoef[u];
					pDataCoef[u] = pDataCoef[v];
					pDataCoef[v] = t;
				}
			}
			t = pDataCoef[k * l + k];
			for (int j = k + 1; j <= l - 1; j++)
			{
				int u = k * l + j;
				if (pDataCoef[u] != 0.0)
				{
					pDataCoef[u] /= t;
				}
			}
			pDataConst[k] /= t;
			for (int j = k + 1; j <= l - 1; j++)
			{
				int u = k * l + j;
				if (pDataCoef[u] == 0.0)
				{
					continue;
				}
				for (int i = 0; i <= l - 1; i++)
				{
					int v = i * l + k;
					if (i != k && pDataCoef[v] != 0.0)
					{
						nIs = i * l + j;
						pDataCoef[nIs] -= pDataCoef[v] * pDataCoef[u];
					}
				}
			}
			for (int i = 0; i <= l - 1; i++)
			{
				int u = i * l + k;
				if (i != k && pDataCoef[u] != 0.0)
				{
					pDataConst[i] -= pDataCoef[u] * pDataConst[k];
				}
			}
		}
		for (int k = l - 1; k >= 0; k--)
		{
			if (k != pnJs[k])
			{
				double t = pDataConst[k];
				pDataConst[k] = pDataConst[pnJs[k]];
				pDataConst[pnJs[k]] = t;
			}
		}
		return true;
	}

	public bool GetRootsetTlvs(Matrix mtxResult)
	{
		int l = mtxLECoef.GetNumColumns();
		mtxResult.Init(l, 1);
		double[] x = mtxResult.GetData();
		double[] pDataConst = mtxLEConst.GetData();
		double[] t = new double[l];
		for (int i = 0; i < l; i++)
		{
			t[i] = mtxLECoef.GetElement(0, i);
		}
		double[] s = new double[l];
		double[] y = new double[l];
		double a = t[0];
		if (a == 0.0)
		{
			return false;
		}
		y[0] = 1.0;
		x[0] = pDataConst[0] / a;
		for (int k = 1; k <= l - 1; k++)
		{
			double beta = 0.0;
			double q = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				beta += y[j] * t[j + 1];
				q += x[j] * t[k - j];
			}
			if (a == 0.0)
			{
				return false;
			}
			double c = (0.0 - beta) / a;
			s[0] = c * y[k - 1];
			y[k] = y[k - 1];
			if (k != 1)
			{
				for (int i = 1; i <= k - 1; i++)
				{
					s[i] = y[i - 1] + c * y[k - i - 1];
				}
			}
			a += c * beta;
			if (a == 0.0)
			{
				return false;
			}
			double h = (pDataConst[k] - q) / a;
			for (int i = 0; i <= k - 1; i++)
			{
				x[i] += h * s[i];
				y[i] = s[i];
			}
			x[k] = h * y[k];
		}
		return true;
	}

	public bool GetRootsetGaussSeidel(Matrix mtxResult, double eps)
	{
		int k = mtxLECoef.GetNumColumns();
		mtxResult.Init(k, 1);
		double[] x = mtxResult.GetData();
		double[] pDataCoef = mtxLECoef.GetData();
		double[] pDataConst = mtxLEConst.GetData();
		double p;
		for (int i = 0; i <= k - 1; i++)
		{
			int u = i * k + i;
			p = 0.0;
			x[i] = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				if (i != j)
				{
					int v = i * k + j;
					p += Math.Abs(pDataCoef[v]);
				}
			}
			if (p >= Math.Abs(pDataCoef[u]))
			{
				return false;
			}
		}
		p = eps + 1.0;
		while (p >= eps)
		{
			p = 0.0;
			for (int i = 0; i <= k - 1; i++)
			{
				double t = x[i];
				double s = 0.0;
				for (int j = 0; j <= k - 1; j++)
				{
					if (j != i)
					{
						s += pDataCoef[i * k + j] * x[j];
					}
				}
				x[i] = (pDataConst[i] - s) / pDataCoef[i * k + i];
				double q = Math.Abs(x[i] - t) / (1.0 + Math.Abs(x[i]));
				if (q > p)
				{
					p = q;
				}
			}
		}
		return true;
	}

	public void GetRootsetGrad(Matrix mtxResult, double eps)
	{
		int k = GetNumUnknowns();
		mtxResult.Init(k, 1);
		double[] x = mtxResult.GetData();
		Matrix mtxP = new Matrix(k, 1);
		double[] p = mtxP.GetData();
		double[] pDataCoef = mtxLECoef.GetData();
		double[] pDataConst = mtxLEConst.GetData();
		double[] r = new double[k];
		for (int i = 0; i <= k - 1; i++)
		{
			x[i] = 0.0;
			p[i] = pDataConst[i];
			r[i] = pDataConst[i];
		}
		for (int i = 0; i <= k - 1; i++)
		{
			Matrix mtxS = mtxLECoef.Multiply(mtxP);
			double[] s = mtxS.GetData();
			double d = 0.0;
			double e = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				d += p[j] * pDataConst[j];
				e += p[j] * s[j];
			}
			double alpha = d / e;
			for (int j = 0; j <= k - 1; j++)
			{
				x[j] += alpha * p[j];
			}
			Matrix mtxQ = mtxLECoef.Multiply(mtxResult);
			double[] q = mtxQ.GetData();
			d = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				r[j] = pDataConst[j] - q[j];
				d += r[j] * s[j];
			}
			double beta = d / e;
			d = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				d += r[j] * r[j];
			}
			d = Math.Sqrt(d);
			if (d < eps)
			{
				break;
			}
			for (int j = 0; j <= k - 1; j++)
			{
				p[j] = r[j] - beta * p[j];
			}
		}
	}

	public bool GetRootsetMqr(Matrix mtxResult, Matrix mtxQ, Matrix mtxR)
	{
		int k = mtxLECoef.GetNumRows();
		int l = mtxLECoef.GetNumColumns();
		if (k < l)
		{
			return false;
		}
		mtxResult.SetValue(mtxLEConst);
		double[] pDataConst = mtxResult.GetData();
		mtxR.SetValue(mtxLECoef);
		double[] pDataCoef = mtxR.GetData();
		if (!mtxR.SplitQR(mtxQ))
		{
			return false;
		}
		double[] c = new double[l];
		double[] q = mtxQ.GetData();
		for (int i = 0; i <= l - 1; i++)
		{
			double d = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				d += q[j * k + i] * pDataConst[j];
			}
			c[i] = d;
		}
		pDataConst[l - 1] = c[l - 1] / pDataCoef[l * l - 1];
		for (int i = l - 2; i >= 0; i--)
		{
			double d = 0.0;
			for (int j = i + 1; j <= l - 1; j++)
			{
				d += pDataCoef[i * l + j] * pDataConst[j];
			}
			pDataConst[i] = (c[i] - d) / pDataCoef[i * l + i];
		}
		return true;
	}

	public bool GetRootsetGinv(Matrix mtxResult, Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps)
	{
		int k = mtxLECoef.GetNumRows();
		int l = mtxLECoef.GetNumColumns();
		mtxResult.Init(l, 1);
		double[] pDataConst = mtxLEConst.GetData();
		double[] x = mtxResult.GetData();
		Matrix mtxA = new Matrix(mtxLECoef);
		if (!mtxA.InvertUV(mtxAP, mtxU, mtxV, eps))
		{
			return false;
		}
		double[] pAPData = mtxAP.GetData();
		for (int i = 0; i <= l - 1; i++)
		{
			x[i] = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				x[i] += pAPData[i * k + j] * pDataConst[j];
			}
		}
		return true;
	}

	public bool GetRootsetMorbid(Matrix mtxResult, int nMaxIt, double eps)
	{
		int k = GetNumUnknowns();
		int i = nMaxIt;
		LEquations leqs = new LEquations(mtxLECoef, mtxLEConst);
		if (!leqs.GetRootsetGauss(mtxResult))
		{
			return false;
		}
		double[] x = mtxResult.GetData();
		double q = 1.0 + eps;
		while (q >= eps)
		{
			if (i == 0)
			{
				return false;
			}
			i--;
			Matrix mtxE = mtxLECoef.Multiply(mtxResult);
			Matrix mtxR = mtxLEConst.Subtract(mtxE);
			leqs = new LEquations(mtxLECoef, mtxR);
			Matrix mtxRR = new Matrix();
			if (!leqs.GetRootsetGauss(mtxRR))
			{
				return false;
			}
			double[] r = mtxRR.GetData();
			q = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				double qq = Math.Abs(r[j]) / (1.0 + Math.Abs(x[j] + r[j]));
				if (qq > q)
				{
					q = qq;
				}
			}
			for (int j = 0; j <= k - 1; j++)
			{
				x[j] += r[j];
			}
		}
		return true;
	}
}
